home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Collection of Tools & Utilities
/
Collection of Tools and Utilities.iso
/
basic
/
chaosexe.zip
/
XINITIAL.TRU
< prev
next >
Wrap
Text File
|
1980-01-01
|
3KB
|
110 lines
!PROGRAM TITLE "XINITIAL( EVOLUTION)"
!THIS PROGRAM DISPLAYS THE EVOLUTION OF A BLOCK OF PHASE POINTS FOR THE PENDULUM
LIBRARY "SGLIB.TRC"
!
DECLARE DEF accel
DIM A(1),B(1)
INPUT prompt"Input driving force strength: ":g
INPUT prompt"Input damping (If no damping then input 9999999):":q
PRINT"THE NEXT SET OF PROMPTS ASK FOR THE BOUNDARIES OF THE AREA OF PHASE"
PRINT"POINTS WHOSE INITIAL EVOLUTION WILL BE GRAPHED."
INPUT PROMPT"LOWEST X VALUE:":LOWX
INPUT PROMPT"HIGHEST X VALUE:":HIGHX
INPUT PROMPT"LOWEST Y VALUE:":LOWY
INPUT PROMPT"HIGHEST Y VALUE:":HIGHY
INPUT Prompt"Input min. and max. time:":tmin,tmax
INPUT prompt"Input TIME STEP AS A FRACTION OF THE FORCING PERIOD: ":phi
LET PHISTEP=PHI*2*PI
CALL PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
CALL SETXSCALE(XMIN,XMAX)
CALL SETYSCALE(YMIN,YMAX)
CALL SETTEXT("EVOLUTION OF INITIAL PENDULUM STATES","THETA","OMEGA")
CALL RESERVELEGEND
!
DATA O,O
CALL DATAGRAPH(A,B,1,O,"RED")
CALL GOTOCANVAS
!
!LOOP FOR BLOCK OF POINT
FOR XINT = LOWX TO HIGHX STEP .1
FOR VINT = LOWY TO HIGHY STEP .1
CALL GRAPHPOINT(XINT,VINT,1)
LET T = 0
LET X=XINT
LET V=VINT
!CALCULATION AND GRAPHING BLOCK
LET phi=phi*2*pi
FOR i=1 to 1000000
CALL rk4(x,v,tstep,xnew,vnew,t,w,g,q) ! Take a step of size tstep
LET tshalf=tstep/2
CALL rk4(x,v,tshalf,xnh,vnh,t,w,g,q) !Take two half steps
CALL rk4(xnh,vnh,tshalf,xn,vn,t+tshalf,w,g,q)
LET d1=abs(xn-xnew)
LET d2=abs(vn-vnew)
LET delta=max(d1,d2)
IF delta<eps then
IF t>tmin then
LET tnew=t+tstep
LET w1=mod(-W*T,PHISTEP) !Check for TIME STEP TO PLOT POINT
LET w2=mod(W*TNEW,PHISTEP)
IF w1<w*tstep then
IF w2<w*tstep then
LET ts=w1/w
CALL rk4(x,v,ts,xp,vp,t,w,g,q) !CALCULATES POINT AT TIME STEP
IF abs(xp)>pi then LET xp=xp-2*pi*abs(xp)/xp
CALL GRAPHPOINT(XP,VP,1)
END IF
END IF
END IF
LET x=xnew
LET v=vnew
LET t=t+tstep !Expand step size
LET tstep=tstep*.95*(eps/delta)^.25
IF abs(x)>pi then !bring theta back into range
LET x=x-2*pi*abs(x)/x
END IF
ELSE !else reduce step size
LET tstep=tstep*.95*(eps/delta)^.2
END IF
IF t>tmax then LET i=1000001
NEXT i
NEXT VINT
NEXT XINT
LET G$=STR$(G)
LET Q$=STR$(Q)
CALL ADDLEGEND("G="&STR$(G)&" Q="&STR$(Q)&" TIME STEP="&STR$(PHISTEP*1.5),0,1,"WHITE")
CALL DRAWLEGEND
END
!
!
SUB rk4(x,v,tstep,xnew,vnew,t,w,g,q)
DECLARE DEF accel
LET xk1=tstep*v
LET vk1=tstep*accel(x,v,t,w,g,q)
LET xk2=tstep*(v+vk1/2)
LET vk2=tstep*accel(x+xk1/2,v+vk1/2,t+tstep/2,w,g,q)
LET xk3=tstep*(v+vk2/2)
LET vk3=tstep*accel(x+xk2/2,v+vk2/2,t+tstep/2,w,g,q)
LET xk4=tstep*(v+vk3)
LET vk4=tstep*accel(x+xk3,v+vk3,t+tstep,w,g,q)
LET vnew=v+(vk1+2*vk2+2*vk3+vk4)/6
LET xnew=x+(xk1+2*xk2+2*xk3+xk4)/6
END SUB
DEF accel(x,v,t,w,g,q)
LET accel=-sin(x)-(1/q)*v+g*cos(w*t)
END def
!
SUB PARAMS(W,EPS,TSTEP,XMIN,XMAX,YMIN,YMAX)
LET W=0.66666666
LET EPS=1.0E-6
LET TSTEP=0.5
LET XMIN=-PI
LET XMAX=PI
LET YMIN=-4
LET YMAX=4
END SUB